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ABSTRACT 

We examine the impact of dark matter particle resolution on the formation of a bary- 
onic core in high resolution adaptive mesh refinement simulations. We test the effect 
that both particle smoothing and particle splitting have on the hydrodynamic proper¬ 
ties of a collapsing halo at high redshift (z > 20). Furthermore, we vary the background 
field intensity, with energy below the Lyman limit (< 13.6 eV), as may be relevant 
for the case of metal-free star formation and super-massive black hole seed formation. 
We find that using particle splitting methods greatly increases our particle resolution 
without introducing any numerical noise and allows us to achieve converged results 
over a wide range of external background fields. Additionally, we find that for lower 
values of the background field a lower dark matter particle mass is required. We define 
the radius of the core as the point at which the enclosed baryonic mass dominates over 
the enclosed dark matter mass. For our simulations this results in R CO re ~ 5 pc. We 
find that in order to produce converged results which are not affected by dark matter 
particles requires that the relationship M core /MDM > 100.0 be satisfied, where M core 
is the enclosed baryon mass within the core and Mdm is the minimum dark matter 
particle mass. This ratio should provide a very useful starting point for conducting 
convergence tests before any production run simulations. We find that dark matter 
particle smoothing is a useful adjunct to already highly resolved simulations. 

Key words: Cosmology: theory - large-scale structure - first stars, methods: numer¬ 
ical 


1 INTRODUCTION 

Gravitational collapse of dark matter has been well 
studied over the past four decades with increas- 


(iDavis et alJ 19851: Frenk et al.l 19881: Warren et al.l 19921: 

Cen et al. 1994: Gelb & Bertschineei 19941: Navarro et al. 

1996 

1997: Hernauist et al. 


19961; Jenkins et al. 

2001 

: Wambsganss et all 2004 


Sm'ingel et all 20081: 

Diemand et al.l 20081: Bosch et al. 

201J) leading to the con- 


firmation that the large-scale structure of the Universe is 
built up through hierarchical collapse of small dark matter 
haloes. It is within the gravitational potential wells of these 
dark matter structures that galaxies and their associated 
active galactic nuclei a n d quasars form ( Whrte_&^ees|ll978|: 
White&_j^renk|_ 1991: Springel et al. 2005|: Croton et alj 


20061: 1 Johansson et al .1120121: IVogelsberger et al.ll2014lh 


Dark matter is modelled as a collisionless fluid and in 
principle can b e solved using the col li sionless Bolt z mann 
equation (e.g . IShlosman et al.l Il979l : iRasio et al.l Il989l : 
iHozumil 1 19971 ). However, the computational demands of 
solving this six dimensional equation means that instead 
dark matter is typically treate d in a Monte-Carlo fa shion 
(however see recent work by lHahn Angulo! 120151 ) and 
represented by a set of collisionless particles which interact 
with each other gravitationally. Solving the gravitational 
interactions of a many-body system is then dealt with in 
one of two ways, direct summation (e.g. von Hoernerlll960l : 
lAarsethlIl963l : Ivon Hoerneil 1 19631 : ISteinmet j|l996l ) ~or more 
commonly, in cosmological simulations, usi ng an approxi ¬ 
mate treat ment involving a Tree algorithm ( Bames_&_JIui 
_198C, Il989|) or a particle mesh alg orithm (lEfstathjou et al 


1A8E: Hockney fe Eastwoodl 1 19881 ) - see also iBertschingerl 


E-mail: j ohn. regan@helsinki. fi 


(1993) f° r a comprehensive review. 

The subsequent addition of gas dynamics to the 
simulations increases the complexity of the simulations sub- 
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2 J.A. Regan et al. 


stantially, requiring not just a detailed treatment of the gas 
(fluid) interaction with the dark matter, but also of a range 
of physical processes involving the gas including radiative 
cooling, star formation and the tracking of individual 
chemical species. Incorporating the gas dynamics into the 
simulation usually also takes the form of one of two mech¬ 
anisms - smoothed particle hydr odynamics (SPH) (iLucvl 


1 19771: Gingold fe Monaghan ll977j) or gr id based methods 
jR^uet^L|_ 1990|: ICen et al Il990l : IStone fe Normanl 1 19921 : 


iBrvan et al 


19951 1 where the hydrodynamical equations are 


solved on the grid using finite difference tech niques possibly 
involving deformin g or adaptive meshes (e.g. ISpringelll20ld : 
IBrvan et alJl2014h . 

Here we focus on the effect the mass resolution of the 
dark matter particles can have on the gas dynamics within a 
collapsing halo where the gas densities can exceed the dark 
matter densities by several orders of magnitude as happens 
in the centres of newly forming galaxies, mini-haloes and 
clusters of galaxies. In order to conduct this case study we 
use the Enzo code (see (|2]for details on Enzo) to follow the 
collapse of gas within collapsing haloes early in the Universe 
at a time when the cooling of the gas is driven by primordial 
processes and in which the gas has not yet been enriched 
with metals from previous episodes of star formation. This 
case is relevant for sim ulations involving the first generation 
of stars (Pop I II stars]j Abel et al]l2002 : Bromm et al .lEfOoij : 


Voshidaet_al. [j^OO^^romm&LarsoJj^OOj^^oshidaet^l 


20061: O’Shea fe Normanl 2008: ^shidaetalj 


200f 


Wise fc Abell 120081: lO’Shea fc Normanl 120081: iTurk et al 


Broinmet^LlJ^OOd; Clnxket^lJJ^Ol^JStacvet^l 


2012|: Susa et all 2014 : Stacv &: Bromua 20141 : 


Hirano et all 120141 1 or for the case of direct collapse 


models of supermassive black holes ( 

Shlosman et al.1 19891: 

Haiman & Loeb 200ll: Oh & Haiman 

2002 

: Bromm & Loebl 

20031: lHaiman 2006: Begelman et al. 20061: 

Wise et al.1 20081: 

Regan & Haehneltl l2009allbl: Regan et al.1 2014al bf) which 


are also thought to form out of pristine gas at high redshift. 

Previous studies of dark matter resolution have indi¬ 
cated the negative numerical effects that dark matter can 
have on th e collapse of baryon s in t he context of galaxy 
formation. ISteinmetz fc Whitel (Il997l l showed that using 
dark matter particles with masses above a critical mass 
limit (their work was focused in the context of galaxy 
formation) introduced spurious two body heating effects 
which they found completely suppressed the effects of 
radiative cooling. Although their numerical tests were 
carried out using the SPH technique they argued that 
the results are valid in the case of grid based approaches 
also given the modelling of the dark matter component is 
similar in both instances. And while other studies have 
noted the effect that dark matter mass resolution p l ays in 
galaxy formation sim ulations (e.g. iMachacek et all l200ll : 
iKaufmann et all I2006T ) no specific study has been made, 
to the best of our knowledge, on the effect of dark matter 
resolution on the collapsing corcQ within an embryonic 


1 The core radius is here defined as the radius at which the bary- 
onic enclosed mass exceeds the dark matter mass for the first time. 
Typically, for our simulations this occurs between 1 and 10 parsecs 
from the point of highest density. This definition is closely related 
to the Jeans length, A j, of the gas for a gas with an isothermal 
temperature of 8000 K and number density, n~lx 10® cm -3 . 


galaxy. This is particularly crucial as correctly accounting 
for the cooling and heating effect becomes extremely 
important in determining the characteristics of the final 
object. In this paper we study this exact effect. 

It is also worth noting that apart from the negative 
impact that poorly resolved dark matter particles may have 
on the baryonic physics an additional numerical artefact 
may occur if a dark matter particle flies into the high 
density region from a lower density region. This “flyby” 
and the perturbing effects of its large mass may then induce 
incorrect physical conditions. While we do not specifically 
deal with this issue in this study it can be alleviated using 
particle smoothing (see <j2~4l) which will necessarily prevent 
adverse effects from a large dark matter mass within the 
baryon dominated core. 

We vary the mass resolution of the dark matter 
particles that surround the collapsing core in an effort 
to identify any spurious heating effects that the dark 
matter has on the collapsing core and identify resolution 
guidelines. Furthermore, we vary the background radiation 
field resulting in the formation (and collapse) of mini-haloes 
and haloes with virial temperatures T > 10 4 K and again 
test for dark matter mass resolution convergence. 

The paper is laid out as follows: in 32] we describe 
the numerical approach used and the various techniques 
employed that can affect the dark matter resolution; 
in m we describe the characteristics of the simulations 
used in this study; in fj4] we describe the results of our 
numerical simulations; in sni we discuss the parameters 
for choosing the correct dark matter particle mass and in 
m we present our conclusions. Throughout this paper we 
assume a st andard ACDM cosmology with the following 
parameters (IPlanck Collaboration et al.ll2014l . based on the 
latest Planck data), 11 a ,0 = 0.6817, fl m ,o = 0.3183, flb.o = 
0.0463, os = 0.8347 and h = 0.6704. We further assume 
a spectral index for the primordial density fluctuations of 
n = 0.9616. 


2 NUMERICAL SETUP 

We have used the publicly available adaptive mesh re¬ 
finement (AMR) code Enz<0. AMR codes are a popular 
tool for studying the gravitational and hydrodynamical col- 
lapse of astrophysical objects (althou gh see for example 
lYoshida et al.ll2006l . Idark et al.ll20llbl & iMaver et al.ll20l3 
among others for collapse simulations conducted using SPH 
techniques) due to their inherent ability to extend over many 
magnitudes in dynamic range. Therefore using Enzo as the 
code with which to perform this study is an obvious choice. 
Furthermore, we use version 3. (0 which is the bleeding edge 
version of the code incorporating a range of new features 
(e.g. Goldbaum et al. in prep). 

Enzo was originally developed by Gr eg Bryan and Mike 
Norman at the University of Illino is (|Brymi_J^jormaij| 
19951. Il997l: iNorman fc Brvanl 1 19991 : lO Shea et al] |2004 
Brvan et al.ll2014h . The gravity solver in Enzo uses an N- 
Body adaptive particle-mesh technique (lEfstathiou et all 

2 http://enzo-project.org/ 

3 Changeset: 6037:lc232317e50d 
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Figure 1. The left hand panel shows the spherically averaged dark matter density profile for a dark matter only simulation of Halo A. 
The output is taken at a redshift of z = 23 and the density is scaled by the mass of the hydrogen atom to enable a direct comparison 
with the gas density. The lines are for a minimum dark matter mass of 122.7 Mq, 9.5 Mq and 0.7 Mq. The right hand panel shows 
the difference between the enclosed mass profile for the same halo at the same output time. The differences are calculated against the 
highest resolution simulation (i.e. the dark matter simulation with M = 0.7 Mq). 


1 19851 : iHocknev fc Eastwoodl Il988l : ICouchmanl Il99ll ) while 
the hydrodynamics are evolved using the piecewise parabolic 
method combined with a non-linear Riemann solver for 
shock capturing. The AMR methodology allows for addi¬ 
tional finer meshes to be laid down as the simulation runs 
to enhance the resolution in a given, user defined, region. 

For our simulations we set the maximum refinement 
level to 18. Refinement is triggered in Enzo when the re¬ 
finement criteria are exceeded. The refinement criteria used 
in this work were based on three physical measurements: 

(1) The dark matter particle over-density, (2) The baryon 
over-density and (3) the Jeans length. The first two criteria 

introduce additional meshes when the over-density ( Ap ) 

' Pmean ' 

of a grid cell with respect to the mean density exceeds 

8.0 for baryons and/or DM. Furthermore, we set the Mini- 

mumMassForRefinementExponent parameter to —0.1 mak¬ 

ing the simulation super-Lagrangian and therefore reducing 
the threshold for refinement as higher densities are reached. 
For the final criteria we set the number of cells per Jeans 
length to be 16 in these runs. 


2.1 Chemical Modelling 

We adopt a nine species model for the chemical network 
used including H, H + , He, He + , He ++ , e _ , H 2 , and H~. 
The gas is allowed to cool radiatively during the collapse. 
Additionally, we incorporate an external radiation field with 
an energy below the hydrogen ionisation edge into the simu¬ 
lations conducted here. Specifically we incorporate radiation 
energies from 0.76 eV up to 13.6 eV which affects the three 
species H~,H 2 ,H(( as described in more detail below. The 
intensity of the background radiation is varied over a range 
of values between 0 and 500 x 10 -21 erg cm -2 s -1 Hz -1 sr -1 
with the effective temperature of the blackbody temperature 
set to T = 50000 K in all cases. We use the variable J 21 as 
shorthand for 10 -21 erg cm~ 2 s -1 Hz -1 sr -1 . We initially 


concentrate our study on radiation fields with 80 J 21 and 
500 J 21 and then in 114.31 we explore the effect of varying 
the external field. The shape of the spectrum is fully de¬ 
scribed by the temperature of the blackbody radiation, T e ff 
and the amplitude of the fluctuation, k J 21 , where k is any 
scalar quantity. We normalise the external radiation field 
at the hydrogen ionisation edge. In our simulations we al¬ 
low the external field to affect the rates of three species - 
H~,H 2 ,Hj~, the rates for the three species are given below 
dAbel et al.lll997l ): 


k 


H- 


47T 

~h 



J21 


5(5, T e ff) cr(E) 
5(13.6, T e// ) E 


kn 2 


= 1.38 x 10 9 J 21 


B(E LW ,T e ff) 
5(13.6, T e ff) 


[s" * 1 ] (1) 


+ 

2 


47T 

~h 



J 21 


B{E, T eff ) a(E) dE 
5(13.6, T e// ) E 


where kx is the photoionisation rate for a given 
species, J 21 is the external field with implied units of 
10 -21 erg cm -2 s -1 Hz -1 sr -1 , cr(E) is the cross section at 
a given energy (frequency) and B(E, T e //) is the blackbody 
spectrum at an energy, E, temperature, T e // and Elw is 
the energy in the Lyman-Werner band (12.8 eV). The in¬ 
tegration limits are in eV and do not exceed the hydrogen 
ionisation threshold of 13.6 eV. 


2.2 Realisations 

Using the initial conditions generator MUSIC dHahn fe Abell 
EoH we generated over 5000 random dark matter reali¬ 
sations with Enzo using a 256 3 * * * * * grid. Each simulation was 
run until a redshift of z = 30 with an inline “Friends of 
Friends” halo finder. Several realisations were identified as 
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Figure 2. Each panel shows a dark matter density projection of the 0.5 kpc (physical) surrounding the collapsing halo at a redshift 
of z = 23. The left hand panel has a minimum dark matter mass of Mom ~ 122.7 M@, the middle panel has a minimum dark matter 
mass of Mom ~ 9.5 Mq, and the right hand panel has a minimum dark matter mass of Mom ~ 0.7 Mg. Each projection is centred 
on the point of maximum density in the simulation with a dark matter particle mass of Mom ~ 122.7 Mg. 


havin g a high sigma peak (a > 4.0 at z > 30, see for ex¬ 
ample Regan et al.l l2014a for details on the definition of a). 
Three realisations were found from the initial dark matter 
runs and a single halo realisation was used to conduct this 
study - labeled Halo A. It should be noted that we have also 
conducted resolution tests on the two other haloes found - 
Halo B and Halo C - however the results were broadly sim¬ 
ilar. We discuss the characteristics of Halo B and Halo C in 
1 14.41 but focus the bulk of this study around Halo A. 

Each realisation was then resimulated with baryons in¬ 
cluded and with different levels of nested grids and hence 
dark matter particle resolution. We ran simulations utilis¬ 
ing one, two and three levels of grid nesting. Each nested 
grid reduces the dark matter particle mass resolution by a 
factor of eight. Memory limitations restrict the number of 
nested grids to a maximum of three. For this study we use a 
box size of 2.0 h -1 Mpc which allows us to find large haloes 
(high sigma peaks) while also restricting the dynamic range 
required by the simulation. This means that the maximum 
dark matter particle resolution that Enzo can achieve in the 
simulation is Mdm ~ 103 Mg. Enzo does not dynamically 
refine the particle mass during the simulation when refining 
its grids. In order to rehne the dark matter particles further 
we employ a Particle Splitting procedure. 


2.3 Particle Splitting 

In order to increase the resolution of the dark matter parti- 
cle s stil l further we emp loy t he pa rticle splitting techniques 
of iKitsionas fc Whitworthl (120021 '). several other authors 
have previously adopted the same technique dYoshida et al.l 
l2006l : lGreif et alll201ll : iKim et ahll201ll : iHirano et al.ll2014h . 

We split particles at a redshift of z = 40 so as to minimise the 
effects of any noise introduced by the splitting technique. In 
particle splitting, the distribution of matter within a region 
assumes that the region has a uniform mass distribution and 
so no new perturbation modes are added. Splitting the par¬ 
ticles at a redshift of 40 ensures that we approximate this 
case and furthermore we test for convergence as described 
below. 


Particles are either split once or twice. Each splitting 
reduces the particle mass by a factor of 13 producing 13 
new child particles in the process and so also increasing 
the computational demands of the simulation. We follow 
iKitsionas fe Whitworthl (120021 1 in choosing to split each par¬ 
ent particle into 13 child particles, the choice of 13 is some¬ 
what arbitrary but attempts to take into account two com¬ 
peting arguments - 1) there should not be too large a dif¬ 
ference between the child mass and the parent mass and 2) 
the collective density distribution of the family of children 
should approximate the spherically symmetric density dis¬ 
tribution of the parent. 

We test for convergence using both one and two itera¬ 
tions of splitting. We split the particles in a subregion sur¬ 
rounding the point of maximum density. The location of the 
point of maximum density is found from previous simula¬ 
tions where particle splitting was not used. This allows for 
the increased mass resolution to be targeted in the region 
surrounding the collapse. We select a region of 43.75 li” 1 
kpc (comoving) within which to do the first iteration of par¬ 
ticle splitting. The second splitting, if performed, is done 
within a region half this size - 21.88 h _ kpc (comoving). 
This allows for maximum dark matter resolution alongside 
maximum grid resolution. 

In order to investigate any systematic effect on the dark 
matter particle properties we initially conducted dark mat¬ 
ter only simulations where we split the particles at a redshift 
of z = 40 and subsequently looked at the highest density re¬ 
gion approximately 80 Myrs later at a redshift of z = 23. 
In the runs with baryons the collapse occurs at a redshift of 
approximately z = 23 or earlier and so examining the form¬ 
ing dark matter halo at a redshift of z = 23 is appropriate. 

In Figure[T]we show the dark matter density profiles and 
enclosed mass profiles for Halo A at a redshift of z = 23. The 
profiles are all centred at the point of maximum density as 
found in the simulation where particle splitting was not con¬ 
ducted (i.e. the simulation with a dark matter particle mass 
of Mdm = 122.7 Mg). In the left hand panel the density 
profile shows clearly that the differences between the run 
with no particle splitting (blue line) and the runs with par¬ 
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Table 1. Halo Details 


Halo Name” 

MTot b 

M D m c 

MBaryon d 

^200 

T f . 

vir 

g 

Utt 

Jti,max 

Halo A 

2.3 x 10 7 

2.0 x 10 7 

3.0 X 10 6 

0.35 

9718 

1.05 X 10 9 

Halo B 

2.3 x 10 7 

2.0 x 10 7 

3.0 X 10 6 

0.34 

9098 

1.01 x 10 9 

Halo C 

1.4 x 10 7 

1.2 x 10 7 

1.6 x 10 6 

0.30 

7402 

1.18 x 10 9 


Notes: The above table contains the simulation name”, the total mass 9 (gas & dark matter) at the virial radius 4 [Mg], 
the enclosed mass in dark matter 1 ” at the virial radius [Mg], the enclosed mass in baryons d at the virial radius [Mg], the 
virial radius e [kpc], the virial temperature-^ [K], and the maximum gas number density 9 in the halo [cm -3 ]. All units are 
physical, unless explicitly stated otherwise. 


tide splitting (cyan and magenta lines) are negligible which 
gives us confidence that the dark matter splitting introduces 
little or no additional perturbing effects to the dark matter 
component. This should not be surprising since the splitting 
occurs early in the halo formation process at a time when the 
non-linear effects of the gravitational collapse are still negli¬ 
gible. To further emphasise the point we show the difference 
between the enclosed mass profiles for the same outputs. We 
normalise against the the highest resolution simulation (i.e. 
the simulation with a dark matter particle mass of Mdm 
= 0.7 Mg). The differences are negligible at large radii (R 
> 100 pc) and increase to only a factor of < 0.3 at small 
scales. 

Figure [2] shows a density projection of the dark matter 
density along the x-axis for this halo (Halo A) for each case. 
The left hand panel has a minimum dark matter particle 
mass of Mdm = 122.7 Mg, the middle panel has a mini¬ 
mum dark matter mass of Mdm = 9.5 Mg and the right 
hand panel has a minimum dark matter mass of Mdm = 0.7 
Mg. The halo becomes smoother as the dark matter reso¬ 
lution is increased but otherwise remains visually the same. 
The spatial scale in each case is given in the left panel and 
corresponds to 250 pc physical. 


2.4 Particle Smoothing 

In regions of high density (resolution) the cell size can be¬ 
come very small and the particles can begin to experience 
spurious two body effects. In order to avoid such spurious 
interactions the dark matter particles can be “smoothed” 
onto the grid once the cell size drops below a given thresh¬ 
old. In Enzo this is controlled via the MaximumParticleRe- 
finementLevel parameter. Once the cell size drops below a 
user defined threshold the dark matter mass is interpolated 
onto the grid with some smoothing length h (which will al¬ 
ways be quoted in comoving units). While computationally 
inefficient to do, it is done in only a very small region and 
so its effect on the overall computational time is small. The 
dark matter now effectively behaves as a continuous, rather 
than discrete, distribution in these very high density regions 
and the effects of (large) dark matter particles coming into 
close contact is alleviated. In this study we set the threshold 
value of h to 2.8 comoving parsecs (equivalent to a refine¬ 
ment level of 12). We run simulations where the smoothing 
is initiated and also where the smoothing is not activated. 
In the runs without the smoothing the smoothing length is 
effectively the minimum cell size (i.e. the cell size at maxi¬ 
mum refinement level) which in the case of our simulations 
is 0.04 comoving parsecs at level 18. 


2.5 External Radiation Fields 


As noted above we concentrate our study on ex¬ 
ternal radiation fields with values of 80 J 21 and 

500 J 21 both of which may be relevant for study¬ 
ing the direct collapse of super-massive black holes 


20061: Diikstra et al. 2008: Reean & Haehnelti 2009b;a: 

Shane et al.1 120101: Johnson et alJ [20131: 

Latif et al. 2013. 

2014alltJ:fAearwal et al.l2014: Reean et al.l 

2014®. A value 


of 80 J 21 sets the field to where H 2 formation is inhibited 
but not completely destroyed, the collapse of the halo there¬ 
fore takes place in an environment where the collapse is not 
isothermal. A value of 500 J 21 on the other hand results 
in an isothermal collapse. By examining both scenarios the 
effects of different dark matter resolutions on the direct col¬ 
lapse model can be fully assessed. 

We further assess the impact of smaller radiation fields 
on the collapse. We use field strengths of both 1 J 21 and 
10 J 21 and we also test the case where no background ex¬ 
ists (i.e. 0 J 21 ). All of these fields may be relevant for the 
case where metal-free stars form (e .g. iBromm et al. I ll999l : 
I Abel et akll2002l ; IBromm et al.ll2002l) . A comparison of the 
affect of different field strengths is discussed in 44.31 

We do not examine the affect of changing the H 2 three 
body fo rmation rates in th is work and instead use the 
rates of iMartin et alJ ll 19961 ) througho ut. This topic has 
bee n examined compreh ensive ly by both Turk et all (l2011al) 
and iBovino et al] (l2014l l . The IMartin et al.l ( 19961) rates are 
found to prov i de the most consistent fit from the work 
of iTurk et ahl (12011a are densi ty dependent and follow 
closely the rates of Forrevl (120131 ) which are advocated by 
IBovino et al.l (120141 ). 


3 SIMULATION DETAILS 

The primary goal of this study is to evaluate the impact, if 
any, of the dark matter resolution on the results of the col¬ 
lapse of the gas. The general halo characteristics of Halo A, 
Halo B and Halo C are given in Table [T] The study focuses 
on Halo A for the most part with the characteristics of Halo 
B and Halo C only discussed in 44.41 

Each (re)simulation of Halo A differs in either the min¬ 
imum dark matter particle size or the minimum smoothing 


4 The virial mass is defined as 200 times the mean density of the 
Universe in this case. 
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Radius (pc) 


Radius (pc) 


Figure 3. Halo A: In the top two panels no smoothing of the dark matter particles is activated and so the smoothing length is the 
same as the minimum cell size - 0.04 pc (comoving). In the bottom two panels the smoothing is activated at a refinement level of 12 
(2.8 pc comoving). The bottom right panel shows the enclosed dark matter mass divided by the enclosed baryon mass - MoM/Mearyon- 
See text for further details on the plot characteristics. 


length of the dark matter particle. Furthermore, the exter¬ 
nal radiation field that a simulation is exposed to is varied 
between 0 J 21 and 500 J 21 . The simulations are automati¬ 
cally halted once the maximum refinement level is reached. 
Therefore each simulated halo, which is exposed to the same 
external radiation field, also has a very similar maximum 
density and comparisons are made at this maximum refine¬ 
ment level output. Each simulation is therefore compared at 
a very similar stage in its evolution. 

The smoothing lengths for each run are set either at the 
maximum refinement level which corresponds to 0.04 (re¬ 
finement level 18) comoving pc or they are smoothed at a 
refinement level of 12 which corresponds to 2.80 comoving pc 
(see 1)2.41 for more details on setting the particle smoothing 
lengths). Each parameter used in each simulation is clearly 
marked in the left hand plot window of each figure. To pro¬ 
mote clarity we refer to a simulation with a dark matter 
particle resolution of 0.6 M@ as S06, a simulation with a 
dark matter particle resolution of 8.0 M@ as S8, etc. 


4 RESULTS 

4.1 Halo A - Varying the Mass Resolution and 
Smoothing Lengths 

We begin by examining the temperature profile for the col¬ 
lapsed haloes for the different dark matter mass resolu¬ 
tions and smoothing lengths. In the following cases the field 
strength is set to either 80 J 21 or 500 J 21 , we will examine 
the case of lower background radiation intensities in d4.3l 
In Figure [3] we show the results for Halo A when the 
external field is fixed at 80 J 21 . The top row in the panel 
shows simulations where no smoothing of the dark mat¬ 
ter particles occurs. The bottom panel shows simulations 
where the smoothing is activated at a refinement level of 12 
(2.8 comoving parsecs). Each simulation was run with differ¬ 
ent minimum dark matter particle resolution ranging from 
Mdm = 6646 Mg to Mdm = 0.6 Mg. Looking at the left 
hand column where we have plotted the temperature pro¬ 
files, centred on the densest point in the simulation, we see 
that down to approximately 1 parsec from the centre the 


© 0000 RAS, MNRAS 000, 000-000 



















7 




CD 

f-H 

cd 

Dh 

CD 

Oh 

a 

<D 

H 




Radius (pc) Radius (pc) 


2 

1 

0 

-1 

-2 

-3 


1 

0 

-1 

-2 

-3 


O 

• i-H 

Ptf 

in 

n 

c3 

T3 

02 

m 

O 

73 

pi 

H 


o 

• i—H 

o3 

Ptf 

in 

in 

o3 

"a 

02 

02 

O 

73 

pi 

H 


Figure 4. Halo A: Same as Figure [3] except the external radiation field is set to 500 J 21 . 


smoothed and non-smoothed simulations agree very well. 
Within the core of the halo differences begin to appear. In 
the non-smoothed simulations (top left panel) the two low¬ 
est resolution simulations show a hotter core compared to 
the higher resolution simulations (i.e. those with a particle 
resolution Mdm < 103 Mg). In the smoothed runs the low¬ 
est resolution simulation has a temperature in the core of 
T core ~ 6000 K while all the other simulations converge to a 
core temperature T cor e ~ 2000 K - note also that S8 shows 
some significant heating within the core at a radius of ap¬ 
proximately 0.05 pc. It should also be noted here that we 
do not include the effects of radiation self-shielding in our 
simulations which may be important in the core of the halo. 

The right hand column yields a partial explanation for 
this. Considering first the non-smoothed runs in the top 
right panel - in both S830 and S6646 the core is dark mat¬ 
ter dominated. The right hand panel plots the ratio of the 
enclosed dark matter mass to the enclosed gas mass. Val¬ 
ues greater than 1.0 indicate regions where the dark matter 
dominates compared to the gas. In simulations S830 and 
S6646 the region inside approximately 1 parsec is dark mat¬ 
ter dominated (in S6646 a single dark matter particle sits 


approximately 1 parsec from the maximum gas density but 
its mass dominates the region). The results from S830 and 
S6646 are therefore likely to be spurious as the dark matter 
resolution in each case is insufficient and the dark matter 
mass ends up dominating the central potential. The heating 
of the core, seen in the left hand panel, should therefore be 
treated with caution due to the presence of perturbing forces 
from these large dark matter particles. 

For the case of the smoothed simulations we again plot 
the dark matter to gas enclosed mass ratios in the bottom 
right panel. In this case we again see that S6646 is dark mat¬ 
ter dominated and that S830 is only marginally gas domi¬ 
nated. The higher resolution simulations (S103, S8 and S06) 
on the other hand are all strongly gas dominated. Again in 
this case S6646 and S830 should be treated with caution 
(even though S830 shows converged results). 

Regarding the temperature dichotomy within the core 
seen in the left hand column, the result is not as bad as it 
first appears and reflects the time at which we choose to ex¬ 
amine the core. Excluding S830 and S6646 we see that S06, 
S8 and S103 all agree quite well (with the possible exception 
of the smoothed S8 simulation which shows a bump at ap- 
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proximately 0.05 pc). However, allowing the simulations to 
run at the maximum refinement level for a short time longer 
results in the core heating up in each case. This results from 
the collisional dissociation of H 2 which pushes the core onto 
the atomic cooling track once again - this is shown in the left 
panel of FigureOand discussed in more detail in Appendix lAl 
(as it lies somewhat outside the scope of the current work). 
We further caution the reader that running at the maximum 
refinement level means we are artifically preventing collapse 
which may induce unphysical results and that further re¬ 
search needs to be done to examine this scenario in more 
detail. 

Strong convergence between S06, S8 and S103 is there¬ 
fore achieved, when the external field is set to 80 J 21 , but the 
study indicates that the lower resolution simulations dissoci¬ 
ate the H 2 faster compared to the higher resolution simula¬ 
tions, driven by spurious dark matter effects, which impact 
the hydrodynamics. 

In Figure [4] we show the results when the background 
radiation is set to 500 J 21 with all other parameters remain¬ 
ing unchanged. In this case the external radiation field is 
strong enough to reduce the H 2 fraction by several orders of 
magnitude within the halo. The left hand column of Figure 
[|] shows the temperature profiles which show highly con¬ 
verged results for both the non-smoothed (top panel) and 
smoothed (bottom panel) for even the lowest dark matter 
resolutions. The right hand column shows the dark mat¬ 
ter to gas enclosed mass ratio. The S6646 simulations are 
dark matter dominated in both cases (and yet the temper¬ 
ature profiles remain converged) but the S830 simulation is 
more gas dominated in the smoothed simulations compared 
to the non-smoothed case. The higher resolution runs also 
show stronger convergence when smoothed. 


4.2 Halo A - Comparing Smoothed and 
non-Smoothed Simulations 


The effect of smoothing the particles onto the grid is most 
noticeable in the low resolution simulations (i.e. S6646 & 
S803) while the effects become less pronounced once the 
dark matter particle mass drops or equivalently once the dis¬ 
cretisation velocity, Vdisc, becomes less than the local sound 
speed (equal to a few km/s at these temperatures). We de¬ 
fine Vdisc roughly as the velocity that a dark matter particle 
can impart onto its surroundings 

v„_- VTf- (2) 

where G is the usual gravitational constant, Mdm is the mass 
of a single dark matter particle and R s is the smoothing 
length. This relationship can also be recast in terms of a 
discretisation temperature resulting in the equation 


Tdisc = 


^niHGMDM 


Rs kb 


(3) 


where fi is the mean molecular weight, set to be 1.22 in 
this study, and kb is the Boltzmann constant. This equation 
can be further expanded by replacing Rg with max(R, Rg) 
where R is any radius becoming 


Tdisc (R) 


/iniHGMDM 
max(R, Rs) kb 


(4) 


Equation [4] can then be used to determine the radius at 
which, for a given dark matter particle mass, the dark mat¬ 
ter (under-)resolution can begin to influence the baryon tem¬ 
perature. We will return to this point in 1)4.31 

If we compare the non-smoothed simulations to the 
smoothed ones when the radiation field is set to 80 J 21 (i.e. 
upper and lower rows of Figure [3]) we see that the tempera¬ 
ture profiles in the left hand panels both show a significant 
scatter. The mass ratios are qualitatively similar (right hand 
panels) in both the smoothed and non-smoothed runs with 
the higher resolution (S103, S8 and S06) well converged but 
with S6646 and S830 showing significant scatter between the 
smoothed and non-smoothed runs. For example, the S830 
run in the smoothed simulation is marginally gas dominated 
while in the non-smoothed case S830 is strongly dark matter 
dominated. The overall effect of smoothing when the radi¬ 
ation field is set to 80 J 21 is therefore marginal, the higher 
resolution simulations are already well converged and so 
smoothing only improves the convergence slightly. The lower 
resolution runs are dark matter dominated and smooth¬ 
ing cannot prevent this although smoothing may somewhat 
lessen the adverse effects. 

In comparing the cases where the external radiation 
field is stronger with J set to 500 J 21 the effects are only 
slightly more encouraging. I 11 the 500 J 21 case the gas is 
warmer and hence the sound speed is correspondingly higher 
meaning that the discreteness effects, quantified by Vdisc are 
somewhat lessened. This has the effect that for the temper¬ 
ature profiles the results appear completely converged as 
the effects of the radiation overwhelm any spurious dark 
matter effects. Nonetheless in both the 80 J 21 case and the 
500 J 21 case the S6646 simulation is strongly dark matter 
dominated and smoothing the dark matter onto the grid at 
this stage only alleviates slightly this particular ill effect. 
However, in the case of S830 we observe that while in the 
non-smoothed case the core is marginally gas dominated it 
is significantly gas dominated in the smoothed case - al¬ 
though the ratio remains above 0.1. Finally similar to the 
80 J 21 case the simulations with a particle mass Mdm < 103 
Mq are relatively well converged in both the smoothed and 
non-smoothed cases although the smoothed case does show 
less deviation at the smallest scales. 

Therefore, the result is that smoothing tends to improve 
convergence only in isolated cases or when the dark matter 
resolution is already high. However, smoothing poorly re¬ 
solved dark matter particles onto the grid will likely have 
no positive effect and it is therefore more important to ad¬ 
equately resolve a given region with a sufficiently high dark 
matter mass resolution and then to smooth the particles 
onto the grid to remove any lingering discrete effects from 
the dark matter particles. 


4.3 Halo A - Changing the External Field 
Amplitude 

So far we have selected only external fields of 80 J 21 or 500 
J 21 . These fields are relevant for investigating the case of 
a collapsing halo exposed to a nearby strong UV source as 
may be required in the direct collapse fo rmation model (e.g. 
iDiikstra et alj 120081 ; IVisbal et al.l l2014f) . We now turn our 
attention to investigating the effects of other external field 
amplitudes which may be expected at this redshift. In par- 
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Figure 5. Halo A: Each row shows the temperature profile and the ratio of the dark matter enclosed mass to the baryon enclosed 
mass. The y-axis label is simply identified as “Ratio” for convenience. The top row shows the case where the particle mass is set to 
Mdm = 103 Mg and no smoothing is employed. The dark matter mass resolution then varies between each row with every other 
parameter remaining unchanged. In each individual panel the external field varies between 0 J 21 and 500 J 21 . The dotted line in the 
left hand panels is the “discretisation temperature” as defined in equation [T] At a given scale it shows the effect a dark matter particle 
of a given mass can have on the baryon temperature. 


ticular, we choose three further field amplitudes of 0 J 21 (i.e. 
no background radiation), 1 J 21 and 10 J 21 . To test for any 
discreteness effects from under-resolved dark matter parti¬ 
cles we do not smooth the particles onto the grid so the 
smoothing length is set to the cell size of the finest grid - 
0.04 comoving parsecs in each case. 

In Figure [5] we have plotted the radial profiles of the 
temperature and the dark matter mass to gas ratio with the 
mass resolution increasing between each panel as indicated 
in the right hand column of the figure. 

Examining first the top panel where the dark mat¬ 
ter particle resolution is set to Mdm = 103 M@ we see 


that the simulations with lower field amplitudes (0 J 21 , 1 
J 21 and 10 J 21 ) all show a strong temperature increase to¬ 
wards the core of the halo. As the mass resolution is in¬ 
creased (middle and bottom rows) this phenomenon disap¬ 
pears and the temperature profiles converge to temperatures 
close to T ~ 1000 K as would be expect ed for the forma¬ 
tion of metal-free stars at this redshift (|Abel_et_aL 2002|; 
Yoshidn L et_aL[ 20061. 20081: Bromm et alj 2009||_ &nith_etjd] 


20091: iTurk et al.l|2009| : IStacv et al.ll2010l;Tciark et al.ll2011al : 

Stacv fc Brommfhoil 


In the lowest mass resolution simulations this spurious 
behaviour can be explained by examining the dark matter 


© 0000 RAS, MNRAS 000, 000-000 












































10 


J.A. Regan et al. 




<D 

f-H 

cd 

Dh 

CD 

Oh 

a 

0) 

H 




10 4 


10 3 


10 2 


10 4 


10 3 


10 2 


10 1 




<u 

f-i 

cd 

f-i 

0) 

a 

a 

0) 

H 




0) 

f-i 

3 


cd 

Jh 

0) 

a 

a 

0) 

H 


Radius (pc) 


Radius (pc) 


Figure 6. Halo B & Halo C : The temperature profile for haloes B & C are shown for two different dark matter particle resolutions. 
Halo B is shown in the top row while halo C is displayed in the bottom row. The particle resolution is set at Mdm = 103 Mg in the 
left hand column while it is set to Mdm = 0.6 Mg in the right hand column. The dotted line in each panel is our analytical estimate 
of the discretisation temperature, Tdisc (H). 


to gas ratio in the right hand column. In the lowest mass 
resolution simulations (Mdm = 103 Mg, top right panel) 
the simulations with external radiation amplitudes of 0 J 21 , 
1 J 21 and 10 J 21 are all dark matter dominated (or only 
mildly gas dominated) within the core of the halo resulting 
in spurious heating of the gas due to the presence of large 
dark matter particles (we saw this before in Figure[3]for ex¬ 
ample). As we increase the dark matter particle resolution 
we see this behaviour disappear as the dark matter to gas 
ratio drops significantly and the temperature profiles con¬ 
verge to T~ 1000 K as expected. 

I 11 the temperature panels (left hand panels) we have 
also overplotted our analytical estimation of the discreti¬ 
sation temperature, Tdisc(R), (see equation [H. This gives 
the scale at which we expect the dark matter particle mass 
to start to impact (negatively) on the temperature. In the 
top left panel we see that the large dark matter particle 
mass used here (Mdm = 103 Mq) means that the Tdi sc (R) 
function intersects the 0 J 21 profile very close to where the 
baryon temperature starts to show a spurious increase. A 


similar effect is seen in the middle left panel where the par¬ 
ticle mass is set to (Mdm = 8 Mq). In the bottom left panel 
the Tdi S c(R) function does not intersect any of the temper¬ 
ature profile lines. This indicates that at this particle res¬ 
olution we would not expect any impact from dark matter 
heating effects at the spatial scales probed in our simula¬ 
tions and this is exactly what we see. 

Taking the field strength of 1 J 21 as the fiducial case we 
find that the core (< 1 pc) of the simulation was resolved by 
only 36 dark matter particles (with a Vdisc ~ 18 km/s and 
Tdisc ~ 47000 K) when the mass resolution was set to Mdm 
= 103 Mq. When the resolution was increased to Mdm = 
8 Mq the number of dark matter particles in the core in¬ 
creased to 473 particles (Vdisc ~ 5 km/s and Tdisc ~ 3600 
K) and increasing the resolution to Mdm = 0.6 Mq resulted 
in the core being resolved by 3072 particles (Vdisc ~ 1 km/s 
and Tdisc ~ 300 K). 

The smaller haloes, those with a gas core less than ~ 10 4 
Mq , which collapse when the external field amplitude is rel¬ 
atively small require a much higher dark matter resolution. 
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Figure 7. Halo A: Radial profiles of the enclosed baryonic mass divided by the mass of the (minimum) dark matter particle mass, 
MdM) used in each simulation. In order to satisfy our resolution requirement the ratio of the enclosed baryonic mass to the dark matter 
particle mass should exceed 100.0 at a radius of 1 parsec. In the left hand panel the external radiation field is set to 1 J 21 while in the 
right hand panel the radius is set to 500 J 21 . The dashed lines marks the radius at 1 pc (vertical dashed line) and the ratio value of 
100 (horizontal dashed line). 


In this case we find that the dark matter particle size must 
be less than ~ 8 Mq (although even at Mdm ~ 8 Mq we 
see some spurious heating of the core for the case of J = 0 
J 21 ) and ideally less than ~ 1 Mq in order to ensure that 
the core of the halo is dominated by baryons and not by 
spuriously large dark matter particles. 


4.4 Comparing Halo B and Halo C 

We now briefly examine the hydrodynamic properties of 
Halo B and Halo C. In Figure [6] we have plotted the tem¬ 
perature profiles of Halo B and Halo C. Halo B and Halo C, 
similar to Halo A, represent high sigma peaks in the density 
field. As a result they contain a halo which collapses early 
even when exposed to a strong external radiation field i.e. 
the collapse redshift when the external field is set to 500 J 21 
is z co i ~ 23.1 for Halo B and z co i ~ 23.2 for Halo C. The top 
row of Figure [6] contains plots from Halo B and the bottom 
row is for Halo C. The mass resolution is set at Mdm = 
103 Mq in the left hand column and Mdm = 0.6 Mq in the 
right hand column as indicated. No smoothing of the dark 
matter particles is employed. The external field is varied in 
each panel to show the effect that the external radiation field 
and the dark matter resolution has on the forming halo in 
each case. Finally, we also overplot Tdi SC (R) (equation Q]) as 
a dotted line in each panel. 

Both haloes display very similar properties to Halo 
A (Figure [ 0 ]). For low field strengths (J < 10 J 21 ) and 
Mdm = 103 Mq the temperature increases strongly within 
the core of the halo resulting in temperatures close to T 
~ 10 4 K. This spurious heating in the core of the halo is 
caused by a dark matter particle mass that is too coarse. 
The large dark matter particle masses induce perturbing 
forces on the gas which heats the core. This conclusion is 
supported by our analytical estimate of the heating effect 
given by Tdisc(R). In left hand panels, where the mass res¬ 


olution is comparatively poor, the temperature profiles for 
the low J 21 values begin to increase in the core at a radius 
very close to where the analytical estimate predicts that heat 
may be imparted from dark matter interactions. 

However, increasing the mass resolution to Mdm = 0.6 
Mq resolves this spurious heating and in this case the tem¬ 
peratures converge to T ~ 1000 K as seen for Halo A in 
Figure [5] Again we see this through out analytical estimate 
of the temperature as well. In the right hand panels Tdisc (R) 
is always below the temperature of the baryons and hence 
has no effect on the thermal physics at the scales probed 
in these simulations. Similar to the results for Halo A it is 
clearly demonstrated that for small external radiation fields 
and hence core masses M core < 10 4 Mq, a dark matter par¬ 
ticle mass of Mdm ~ 1 Mq is required. 

For field strengths with J > few x 10 J 21 a very high 
resolution dark matter particle is not required and in these 
cases a particle mass of Mdm ~ 100 Mq is adequate due to 
the more massive gas core in these cases. 


5 CHOOSING THE CORRECT DARK 
MATTER PARTICLE SIZE 

Using a suite of simulations we attempted to constrain the 
dark matter resolution required to achieve convergence in 
the baryonic properties of haloes collapsing at high redshift. 
The particle mass required depends only very weakly on the 
halo (we see no differences between the three haloes simu¬ 
lated) but depends strongly on the enclosed baryonic mass 
within the core of the halo (here defined as the radius at 
which the baryonic enclosed mass exceeds the dark matter 
enclosed mass - typically occuring between 1 parsec and 10 
parsecs in our simulations). We used different external radi¬ 
ation fields which has the effect of varying the mass of the 
core that forms and the cooling processes available in the 
collapsing halo. 
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Enclosed Mass Ratio 



Figure 8. Halo A: Radial profiles of the enclosed baryonic mass divided by the mass of the (minimum) dark matter particle mass, 
M d m, used in each simulation. The scaling ratio is this time plotted against the enclosed mass ratio (i.e. enclosed dark matter mass 
divided by enclosed baryon mass). The dashed lines marks the point of baryon - dark matter enclosed mass equality (vertical dashed 
line) and the ratio value of 100 (horizontal dashed line). 


In Figure [7] we have plotted the enclosed baryon mass 
scaled by the minimum dark matter particle mass used in 
each simulation against the radius. We show the results for 
two representative cases; J = 1 J 21 in the left hand panel 
and J = 500 J 21 in the right hand panel. From our previ¬ 
ous results we know that when the external field is set to 1 
J 21 then a dark matter resolution of Mdm < 103.0 Mq is 
required. On the other hand when the field strength is set 
to 500 J 21 then a dark matter resolution of Mdm < 830.0 
Mq is sufficient. 

In the case where J = 1 J 21 we see that the required 
mass resolutions (S8 and S06) have a ratio of the enclosed 
mass within the core, M core , to the dark matter particle 
mass, Mdm, of > 100 at the core radius i.e. M cor e/MDM > 
100.0. Similarly, in the case where J = 500 J 21 we see that 
three dark matter particle masses satisfy the criterion that 
M cor e/MDM > 100.0 and we have found previously that 
these three particle masses (Mdm = 103 Mq, Mdm = 8 
Mq and Mdm = 0.6 Mq) produce converged results for this 
field strength (or equivalently baryonic core mass since the 
accretion rate depends on the core temperature which de¬ 
pends on the radiation field strength). 

In Figure [H] we have plotted the enclosed baryon mass 
scaled by the minimum dark matter particle mass against 
the enclosed mass ratio (i.e. enclosed dark matter mass di¬ 
vided by enclosed baryon mass). The representation clarifies 
that the scaling ratio must exceed a value of 100.0 at the 
point of baryon - dark matter enclosed mass equality (and 
the core must also remain baryon dominated) in order to 
fulfill our resolution requirements as found previously. 

We therefore present this empirically derived result as 
a “rule of thumb” - 

M cor e/MDM > 100.0 (5) 

The minimum dark matter particle mass used must be a fac¬ 
tor of one hundred less than the baryonic mass contained 
within a sphere in which the baryon mass dominates. 


However, we caution that this relationship does not provide 
a sufficient condition. It should instead be used as a start¬ 
ing point for conducting convergence tests for high redshift 
collapse simulations. Our simulations indicate that failing 
to adhere to this relationship will most likely produce spu¬ 
rious results and at the very least this relationship should 
be satisfied and convergence checked before performing full 
production runs. 


6 CONCLUSIONS 

We have conducted a comprehensive study of the dark mat¬ 
ter particle resolution requirements necessary to study the 
collapse of gas at high redshift in the presence of a back¬ 
ground radiation field below the Lyman limit. We have used 
the publicly available Enzo code to conduct the study. Enzo, 
like most other grid (and particle) based codes, represents 
dark matter as discrete particles while solving the hydrody¬ 
namic equations on a Eulerian grid. In this light, we tested 
a number of different parameters that can potentially affect 
the interaction of the dark matter particles with the collaps¬ 
ing baryonic core. 

We investigated the effect of “smoothing” the dark mat¬ 
ter particles onto the grid in regions of very high density in 
order to negate any two body scattering effects that may 
occur when the cell size becomes very small. We probed the 
effect of increasing the dark matter resolution around re¬ 
gions of interest by employing a particle splitting algorithm 
and tested its effectiveness. Furthermore, we varied the back¬ 
ground radiation field from having no external radiation, J 
= 0 J 21 , up to J = 500 J 21 where J 21 is in the usual units 
of lx 10~ 21 erg cm -2 s _1 Hz^ 1 sr -1 . 

We find that smoothing particles onto the grid provides 
some relief from the effects of spurious two body dark mat¬ 
ter interactions and is most effective when the dark matter 
resolution is already comparatively high. With a dark mat¬ 
ter particle mass of Mdm )$ 100 Mq, smoothing the particles 
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onto the grid produces more converged results, when com¬ 
pared against the non-smoothed case, as the dark matter res¬ 
olution is increased but nonetheless having a high dark mat¬ 
ter resolution was the crucial and overriding factor. At more 
intermediate particle masses, Mdm > 100 Mg, smoothing 
can alleviate some of the effects of two-body scattering (com¬ 
pared to the unsmoothed case) but it cannot negate them 
entirely and is somewhat case (halo) dependent. When the 
dark matter mass resolution is low ( i.e. Mdm > 1000 M@) 
smoothing provided no extra positive effect to alleviate spu¬ 
rious effects as expected. 

In summary, smoothing is most effective when the dark 
matter resolution is already high as it wipes out any lin¬ 
gering discreteness effects. It is significantly less effective as 
the mass resolution drops and it should not be seen as a 
fix for under resolved dark matter simulations. As a result 
we advocate the use of smoothing in conjunction with high 
particle resolution only. 

In examining the effectiveness of particle splitting we 
find very encouraging results. We initially tested the effect 
of splitting on the properties of dark matter in dark mat¬ 
ter only simulations and found that the splitting added lit¬ 
tle or no additional numerical noise to the results finding 
that dark matter density profiles and enclosed mass profiles 
remained highly converged both before and after splitting. 
After adding in baryons, splitting the dark matter enabled 
us to remove spurious dark matter effects from otherwise 
under-resolved simulations. 

Using an external radiation field with an energy below 
the Lyman limit we found that the particle mass required 
varied with the strength of the radiation field. This is not 
surprising since a strong radiation field effectively dissoci¬ 
ated H 2 (and its pathway element H“) reducing the cooling 
ability of the halo thus increasing the accretion rate making 
it more massive. We found that as the halo grew the dark 
matter particle mass required also grew (i.e. less dark mat¬ 
ter resolution was required). This means that when one is 
examining the affect that a very strong radiation field has 
on a collapsing halo the dark matter particle mass may be 
greater (by several factors), without any loss in numerical 
accuracy, compared to case where haloes are collapsing in 
the absence of an external radiation field or in the presence 
of a weak background field. 

We find that a useful “rule of thumb” is that the ra¬ 
tio of the enclosed baryonic mass within the core of the 
halo to the individual dark matter mass must exceed 100 
i.e. M COTe /MoM > 100.0. We found that, for the simulations 
conducted here, this relationship must always be satisfied 
to produce converged results. We therefore advise that this 
relationship always be fulfilled and that appropriate con¬ 
vergence tests be conducted in parallel to ensure reliable 
results. 

While this study was carried out using the grid code 
Enzo we expect the results to hold for other grid and non¬ 
grid codes (e.g. SPH codes) that use the N-body technique 
to represent the dark matter component. O n a final note it is 
appro priate to mention the recent work of lHahn fe Angul3 
(1201511 who describe a method to follow the dark matter evo¬ 
lution using a six dimensional phase approach rather than 
the traditional N-body model studied here. Their approach 
dispenses with the discrete nature of the N-body method 
and may help to significantly suppress any discreteness ef¬ 


fects that can arise in modelling collisionless systems in the 
future. 
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A RUNNING AT THE MAXIMUM 

REFINEMENT LEVEL 

As discussed in 44.11 Figure [5] shows the temperature profile 
for Halo A for both a non-smoothed (i.e. li = 0.04 pc, left 
panel) and a smoothed run (i.e. h = 2.8 pc, right panel) 
further into its evolution than presented in Figure [3] The 
simulations are run at the maximum refinement level to ex¬ 
amine the further evolution of the core. As mentioned pre¬ 
viously, running at the maximum refinement level artifically 
prevents collapse while the Jeans length is also not resolved. 
This may potentially induce unphysical results, however we 
are careful, as detailed below, to only run the simulation for 
a few timesteps. In this case we plot only the results for the 
runs in which the collapsing core is gas dominated and there¬ 
fore exclude simulations S6646 and S830 from both plots so 
as to exclude any potential dark matter contamination ef¬ 
fects which may affect the physical results. 

The simulations have been restarted from the last out¬ 
put and allowed to run at the maximum refinement level 
for a short time (~ 130,000 years). No artificial pressure 
floor is employed here so as to exclude any artificial heat¬ 
ing effects although in principle running at the maximum 
refinement level without any ar tihcial pressure suppor t can 
lead to numerical instabilities (jMachacek et all 120011 '). We 
are therefore careful to run the simulations only for a few 
timesteps at this point. As can clearly be seen, as the sim¬ 
ulation is allowed to progress, the temperature in the core 
increases and the core heats up to close to T ~ 6000 K in 
all cases. The reason for this heating is that dense clumps 
of neutral hydrogen form in the core promoting (collisional) 
dissociation of H 2 through the reaction 

H 2 + H —► 3H (6) 

We have used the formulation given in lMartin et al.l (|l996l 'l 
to model this reaction. The dissociation rate then exceeds 
the H 2 formation rate in regions within the core of the halo 
but conversely the formation rate exceeds the dissociating 
rate in other neighbouring parts resulting in the formation of 
hot and cold spots within the core of the halo. The formation 
rate is controlled, at these densities, by the reaction 

H~ + H -> H 2 + e (7) 

In Figure [10] we have plotted four projection plots illus¬ 
trating this effect. The projections are made approximately 
130,000 years after reaching the maximum refinement level 
using a S103 (with no-smoothing) output for illustration (al¬ 
though S8 or S06 could also have been used). In the bottom 
left panel we show the neutral hydrogen number density 
showing the formation of a dense core, in the top left panel 
we show the temperature projection clearly displaying re¬ 
gions of hot and cold gas within the core. The H 2 fraction 
is shown in the top right panel, reflecting the temperature 
projection, we see regions of low H 2 fraction surrounded by 


regions of high H 2 fraction which dramatically affects the 
cooling in these regions and causes bumps in the temper¬ 
ature profile as shown in Figure [9] Finally, in the bottom 
right panel we show the ratio of the H 2 formation rate to 
the H 2 dissociation rate. Again, we see regions of formation 
and dissociation within the core consistent with the other 
projections. The core of the halo is essentially sitting on a 
knife edge between H 2 formation and dissociation. 

It should also be noted that our simulations neglect the 
potentially important effects of self-shielding. However, a 
study of the complete effects of self-shielding on the core 
hydrodynamics is beyond the scope of this paper. 
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Figure 9. Halo A: The left hand panel shows the spherically averaged temperature profile in the centre of the collapsing halo a short 
time 130,000 years) after running at the maximum refinement level. No smoothing of the dark matter particles occurs in the left 
hand panel. The right hand panel shows the same temperature profile for simulations in which smoothing of the dark matter particles 
occurs at 2.8 comoving pc (refinement level 12). 
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Figure 10. Halo A: The four panel projection shows the state of the gas within the core of the halo approximately 130,000 years 
after reaching the maximum refinement level. Each projection is aligned perpendicular to the angular momentum vector. In the top 
left panel is shown the temperature projection showing the central “hot” core surrounded by cooler gas. The H 2 fraction is shown in 
the top right panel and shows that the core contains clumps of molecular gas surrounded by gas where the H 2 fraction is significantly 
lower. The lower right panel shows the ratio of the H 2 formation rate to the H 2 dissociation rate. Green areas represent areas where 
H 2 is forming while blue areas represent areas where the H 2 is being dissociated. Finally, in the lower left panel the H number density 
is shown, in the core of the halo the H number density is at its maximum. 
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